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Abstract. A new method to numerically calculate the nth moment of the spin 
overlap of the two-dimensional ± J Ising model is developed using the identity derived 
by one of the authors (HK) several years ago. By using the method, the nth moment 
of the spin overlap can be calculated as a simple average of the nth moment of 
the total spins with a modified bond probability distribution. The values of the 
Binder parameter etc have been extensively calculated with the linear size, L, up 
to -L = 23. The accuracy of the calculations in the present method is similar to 
that in the conventional transfer matrix method with about 10^ bond samples. The 
simple scaling plots of the Binder parameter and the spin-glass susceptibility indicate 
the existence of a finite-temperature spin-glass phase transition. We find, however, 
that the estimation of is strongly affected by the corrections to scaling within the 
present data (L < 23). Thus, there still remains the possibility that Tc = 0, contrary 
to the recent results which suggest the existence of a finite-temperature spin-glass 
phase transition. 



PACS numbers: 75.50.Lk,02.70.Lq,64.60.Cn,05.50.-Hq 
1. Introduction 

Over the last two decades, investigations of spin-glass problems have been extensively 
performed [1-20]. It is now widely believed that the three-dimensional ±J Ising 
model shows a finite-temperature spin-glass phase transition [1-6], while the critical 
temperature of the two-dimensional ±J Ising model is zero [5-11]. Most of these 
studies have been done using Monte Carlo simulations, where the thermal relaxation 
time in the simulations becomes very large in the low-temperature region. This makes 
the investigations of the two-dimensional models rather difficult, since the calculations 
of the physical quantities have to be performed at very low temperature. In previous 
studies, the data in the finite-size scaling analysis were in good agreement with a 
scaling function with the critical temperature, Tc = 0. The results, however, have not 
completely excluded the possibility of a finite-temperature spin-glass phase transition. 
Recently, Shirakura et a/[12-15] have deduced the existence of a finite-temperature 
spin-glass phase transition of the two-dimensional models using extensive Monte Carlo 
simulations. To clarify the critical properties of the two-dimensional ±J Ising model, 
more precise results in the low-temperature region are necessary. 

The transfer matrix method is free from the problem of thermal equilibration, 
and has been very successfully used to determine the ferromagnetic-nonferromagnetic 
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phase boundary of the two-dimensional ± J Ising model in the p — T plane |l^ 
{p is the concentration of the ferromagnetic bond). But, for the problem of the 
spin-glass phase transition, the transfer matrix method has only been used for the 
calculations of defect energies and correlation functions and has not been 

widely used for direct calculations of the nth moment of the spin ovelap (the spin- 
glass susceptibility, the Binder parameter, etc) since, when we use real replicas in the 
calculations, the maximum linear size applicable is one-half of that in the calculations 
of the nth moment of the total spins. Thus, so far, no extensive result for the spin- 
glass phase transition with direct calculation of the spin-glass susceptibility etc has 
been given by the transfer matrix method. 

In this paper, we present a new method to numerically calculate the nth moment 
of the spin overlap of the two-dimensional ± J Ising model, using the identity derived 
by one of the present authors several years ago jl^ . In this identity, the nth moment of 
the spin overlap is transformed as a simple average of the nth moment of the total spins 
with a modified bond probability distribution. Following a newly developed process, 
explained in section 2, we successively make the bond configurations according to the 
modified bond probability distribution using the Monte Carlo technique. In each bond 
configuration, we calculate the nth moment of the total spins by the transfer matrix 
method. 

We have performed extensive calculations of the nth moment of the spin overlap 
up to the linear size, L — 23. The accuracy of the calculations in the present 
method is similar to that in the conventional transfer matrix method with about 
10^ bond samples. Thus, the statistical errors in the present study are about an 
order of magnitude smaller than those in previous studies. Therefore, we can analyse 
the obtained data in detail using finite-size scaling analysis including the corrections 
to scaling. Our results show that the estimation of Tc is strongly affected by the 
corrections to scaling in the two-dimensional ± J Ising model. Thus, there still remains 
the possibility that Tc — 0, contrary to the recent results by Shirakura a/ |l2, 13, 1^ 
which suggest the existence of a finite-temperature spin-glass phase transition. 



2. New calculation method for the spin-glass order parameter 

We consider the two-dimensional ± J Ising model on an L x i square lattice with only 
nearest neighbour interactions. The Hamiltonian is written as follows: 

n^^Y.^.,s.s„ (1) 

(y) 

where Si = ±1, and the summation of (ij) runs over all the nearest neighbours. We 
take the skew boundary condition in one direction, and the free boundary condition 
in the other direction. Each Tij is determined according to the following probability 
distribution: 

P{nj) = pS{n, - 1) + (1 - p)5{n, + i). (2) 

In this paper, we make J — I. We define the spin overlap, Q, between two replicas in 
each bond configuration as 



N 



(3) 
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where a and /? denote the two rephcas, and N is the total number of spins. When we 
define Kp as 

exp(2Xp) = (4) 
the 2nth moment of the spin overlap is written as 

[< >T^^s^^s^y]p = ^2coshiKp))^^ exp(ifp^Ty) < Q^n >^ ^s^ s^^, (5) 

where < • • • >T,{5°.sf} denotes the thermal average both for the "5"" and " S^" spins 
in a bond configuration, {r} at temperature, T. [■ ■ ■]p denotes the configurational 
average at the ferromagnetic bond concentration, p, and A'b is the number of bonds 
p8[ . By the use of a local gauge transformation, an identity has been derived Q : 

[< >T,s^,.,]. = (2cosh(i.,))^B E -P(i^E-.)§^ < >-.^>-(6) 
where M denotes the total spins, 

N 

M = ys,, (7) 



=1 



< ■ ■ ■ >T.{s} denotes the thermal average for the "S"' spins at temperature, T, and 
Z{K,{t}) is the partition function at the inverse temperature, K{= 1/T), with the 
bond configuration, {r} . We show the summary of the derivation of equation (6) 
in the appendix. (For the details of the derivation, see [19].) When we define the 
modified probability distribution, P2{K, Kp, {t}), for the bond configuration, {r}, as 

we can then write 

[< Q'" >T..{S^..S^}]p = {< >t,{S}}k.K„ (9) 

where {• • -yK.Kp denotes the configurational average by the modified bond probability 
distribution. That is, [< Q^" >t,{s°,S'3}]p temperature, T, with the ferromagnetic 
bond concentration, p, is transformed into the configurational average of < 
M^" >T,{s} by the modified bond probability disribution, P2{K, Kp, {t}). Similarly, 
we can get the following identity : 

[< >T.is}]p = {< A/'" >t,,{S}}k,k,, (10) 

where Tp — 1/Kp . (Note that the above argument applies to any dimension.) 

Hereafter, we explain a new approach to numerically calculate the values of 
[< Q^" >T,{S'=',Si^}]p7 using equation (9). To realize the bond configuration with 
the modified bond probability distribution, P2{K, Kp, {t}), we use the conventional 
Monte Carlo technique. We define W{Tij —Tij , {r}') as the transition probability by 
which the value of the bond, , changes. To guarantee that the stationary probability 
distribution becomes P2{K, Kp, {nj}), the following detailed balance must be satisfied: 

P2{K,Kp,t,,,{t]')W{t,, -n„{T}') 

= P2{K, Kp, -n„ {T]')W{-n, ^ n„{r}') (n) 
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Using equation (8), we obtain 

W{n,^-T,,,{T}') ^ ^ ^2^^ cosh(2JCp) - sinh(2i^p) < n,S,Sj >r„{s} 

Namely, when we can calculate < SiSj >t,{S} in a particular bond configuration, we 
can estimate the transition probability. 

The processes to calculate [< Q^" >t,{s°,S''}]p are as follows: 

1) We start from an arbitrary bond configuration. 

2) Using the transfer matrix method, we exactly calculate the value of < SiSj >t.{s}i 
and successively flip the bond, r^ , according to the transition probability, W{Tij — > 
— r.y , {r}'), using the conventional Monte Carlo technique. 

3) We continue process 2) until the modified bond probability distribution, 
P2{K,Kp,{t}) is realized. 

4) We calculate the value of < M^" >t,{S] for each bond configuration using the 
transfer matrix method. 

5) We repeat processes, 2) and 4). 

6) Finally, the simple average of < Af^" >t,{s} gives the value of [< Q^" >t,{s°,S'3}]p 
with the bond probability distribution, P{Tij). 

We now show the efficiency of this method. We define ria, Ub and ric as the number 
of initial Monte Carlo skip steps, the number of Monte Carlo steps where we calculate 
< M^" >T,{S}^ and the number of independent Monte Carlo runs, respectively. In all 
the calculations, we have evaluated {Tij{0)Tij{t)}K,K.i, using the statistical dependence 
time method and find that the relaxation time of {Tij{Q)Tij{t)} K,Kp is always 
very small even when compared with one Monte Carlo step time. That is, only 
several tens of initial skip steps are enough to realize the stationary bond probability 
distribution. For example, we have compared the values of the spin-glass susceptibility 
Xsg(= [< >T,{s=,S'3}]p/A^) calculated by the present method and that by the 
conventional transfer matrix method using real replicas. Table 1 shows the results 
at T = 0.1 for L = 7. The calculations by the conventional transfer matrix method 
have been done with 10^ independent bond configurations. The error bars of the 
present methods have been estimated in the same way as those of conventional Monte 
Carlo simulations. From table 1, we can see that all the data are consistent within 
the error bars, and the size of the error bars of all the calculations are of the same 
magnitude. Consequently, we find that only 20 steps are enough for the initial Monte 
Carlo skip steps. Furthermore, we have examined whether equation (9) holds or not 
at p = 0.5 — 0.95, T = 0.1 — 0.5 for L = 7. We have also examined whether equation 
(10) holds or not at p = 0.8 - 0.9, T = 0.1 - 0.4 for L = 15. AU the results are 
consistent in a statistical sense, from which we conclude that the present method is 
usable and not affected by systematic errors. 

3. The spin-glass phase transition of the two-dimensional ± J Ising model 

We have extensively investigated the two-dimensional ± J Ising model for p = 0.5. The 
results of the asymmetric case {p > 0.5) will be given in a subsequent paper |2^. For 
p — 0.5, we have calculated [< Q^" >t,{s°,s'='}]p at T = 0.1 — 0.5 with the linear size 
L = 7 — 23. The calculations have been performed with (ua, rii,, ric) = (200, 1000, 100) 
for L<21 and (200,200,480) for L = 23. 
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Table 1. The values of xsG with various (na,ni,,nc) at T = 0.1 for L = 7. The 
value of * is calculated by the conventional transfer matrix method using real replicas 
with 10® bond samples. 



(na,nt,nc) XSG 

(2000,10000,10) 29.076(24) 

(200,1000,100) 29.084(38) 

(20,100,1000) 29.094(35) 

(20,20,5000) 29.045(34) 

* 29.037(35) 



The energy gap between the ground state and the first excitation state is two in 
the unit of the interaction strength. Thus, in finite systems, any physical quantity at 
a very low temperature must saturate to its value at T = 0. We show the temperature 
dependence of the spin-glass susceptibility, xsG) in figure 1. We find, indeed, that the 
values of xsG for each L show the strong saturation near T = 0, and the tendency 
becomes clearer as the system size becomes smaller, as has already been pointed out 
by several authors ||, [T2| ^ . The Binder parameter |Q 

9L = :r(3- ^) (13) 

is widely used for the determination of the critical temperature. The value of the 
Binder parameter becomes asymptotically size independent for large systems at the 
critical temperature. Therefore, the point where this quantity becomes asymptotically 
size independent gives an estimation of the critical point. The simple plot of the 
Binder parameter versus temperature is shown in figure 2. We can clearly see that 
the data for different sizes intersect at almost the same temperature, T = 0.3, and 
the size dependence of the intersection points is very small. We cannot, however, 
immediately conclude that the spin-glass phase transition occurs at T ~ 0.3, since the 
intersection might be due to the strong saturation of the data near T = . 
Therefore, we perform scaling analyses for and XSG- There is no general rule to 
avoid the disturbance from the saturation mentioned above in the scaling analyses. 
Here, we adopt a criterion that every data point is not used all through the scaling 
analyses, when the value of xsG increases less than 3% in the temperature interval, 
AT = 0.05. Although the criterion, which we have determined from the observation in 
figure 1, seems rather artificial, we believe that this criterion systematically removes 
the saturation to T = in a certain sense. Consequently, we use, for example, the 
data with T > 0.35 for L = 7, and with T > 0.2 for L = 23. 

First, we perform the scaling analyses without including the corrections to scaling. 
In this case, the Binder parameter, g^, has the scaling form 

gi = g(Li/^(T-Tc)), (14) 

and that of the spin-glass susceptibility, xsGi is 

XSG = i'-''x(i'/'^(T-Tc)), (15) 

where v is the critical exponent of the spin-glass correlation length, and rj is the critical 
exponent which describes the decay of the correlation at the critical temperature. 
Figure 3 shows the best-fit scaling plot of the Binder parameter, which indicates that 
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Tc ~ 0.3 and the critical exponent v ~ 1.3. We can see that the data at T < Tc 0.3 
and T > Tc fit rather well on one scaling function, which indicates that the spin- 
glass phase transition of the two-dimensional ± J Ising model is a conventional phase 
transition, and there exists a finite long range order at T < Tc. The best- fit scaling 
plot of the spin-glass susceptibility is also shown in figure 4, which indicates that the 
critical exponent t] ~ 0.225. The estimated values of Tc and the critical exponents are 
similar to those determined by Shirakura et a^^. Figures 5 and 6 show the scaling 
plots of the Binder parameter and the spin-glass susceptibility, assuming = 0, 
V — 2.% and r\ = 0.2 where we clearly see the systematic deviations. Namely, 
the scaling plots without including the corrections to scaling strongly indicate the 
existence of a finite-temperature spin-glass phase transition. 

Next, we perform the scaling analyses including the corrections to scaling. We 
take the scaling forms of the Binder paremeter and the spin-glass susceptibility as 

ffL = ff(Li/^(T-Tc))(l + -^), (16) 

and 

XSG = i^-''x(i'/'^(T - Tc))(l + -^). (17) 

There is little quantitative change in figures 3 and 4, even though we fit the data using 
equations (16) and (17). Thus, when we assume Tc = 0.3, the effect of the corrections 
to scaling is rather small. On the other hand, assuming that Tc = 0, // = 2.6 and 
?7 = 0.17 the data of the Binder parameter and the spin-glass susceptibility fit very 
well on one scaling function, respectively, as shown in figures 7 and 8 with lo = lo' = 0.5 
and a = b = —0.3, although the data with a small linear size, L — 7, deviate from the 
scaling function. (To fit the data, we use rj = 0.17, which is, for example, consistent 
with the result in [8], 77 = 0.2 ± 0.05.) Thus, including the corrections to scaling, both 
Tc = and Tc ~ 0.3 are consistent with the scaling analyses. Furthermore, we find 
that every temperature for < T < 0.3 might become the critical temperature, Tc, 
in this scaling form. Consequently, we find that the estimation of the value of Tc is 
strongly affected by the corrections to scaling in the two-dimensional ± J Ising model 
within the present data {L < 23). 

4. Conclusions 

We have developed a new method to numerically calculate [< Q^" >T,{s°,S'3}]p of 
the two-dimensional ±J Ising model, where, using a local gauge transformation, 
[< Q^" >T,{s°,S'9}]p can be calculated as a simple average of < M^" >t,{s} with 
a modified bond probability distribution. By using the present method, we have 
extensively calculated the values of [< Q^" >t,{s°,si^}]pj where the statistical errors 
become about an order of magnitude smaller than in previous studies, and we have 
investigated the scaling analyses including the corrections to scaling. By using the 
scaling analyses without including the corrections to scaling, our data strongly indicate 
a finite-temperature spin-glass phase transition. We find, however, that the estimation 
of Tc is strongly affected by the corrections to scaling within the data with L < 23, 
and that every temperature for < T < 0.3 might be able to become the critical 
temperature. Consequently, our results indicate that there still remains the possibility 
that Tc = 0, contrary to the recent results of Shirakura et al ||l^, |l^, |l^ which suggest 
the existence of a finite-temperature spin-glass phase transition. 
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Appendix A. 

In this appendix, wc briefly explain the derivation of equation (6). We show that both 
sides of equation (6) coincide with each other. 

Using equation (5), [< Q^" >t,{S'^,S0}]p is written as 

1 ^ 
[< Q'" >T,{s»,s.}]p = ^ E exp(if,5]r,,)<(E5r5ff">T,{s»,sn, (Al) 

TiJ = ±l (ij) i=l 

and we abbreviate (2cosh(iirp))^^ as C from now on. 

Here, we perform the following local gauge transformation: 

Tij TijdiGj, Sf Sfdi, Sfai, {ai = ±1) (A2) 

where each ai arbitrarily takes -|-1 or —1. Since < {J2f=i SfS^Y^ >t,{S",SP} is 
invariant under this transformation, we obtain 

1 ^ 

[< >T,{5o,S'3}]p =c ^MK^y^iiCJiaj) < iJ^S^Sff^ >T,{S'^,s^}, (A3) 

r.,=±l (ij) ,;=i 

where we note that the summation over Tijaiaj = ±1 is eqivalent to that over Tij = ±1. 
As each u, arbitarily takes -|-1 or —1, therefore, we take all the summations of u, and 
divide by 2^, namely 

1 1 ^ 

[< Q'" >T,{S-,S'3}]p = Yl 12 ''MKpY^TijaiCTj) < {Y^S^Sff"" >T,{S'',S0} 

o'j=±lTij=±l (ij) i=l 

~C ^ 2^ <(2^'^i'^i) >T,{s°,s/3} (A4) 

Tij=±l i=l 

Next, we consider the r.h.s. (we denote it as A) of equation (6). The r.h.s. of 
equation (6) is written as 



rij=±l (ij) V ' I J/ i=i 

- ' V c::p(/- V r ) ^">) ^^^^^ exp(j^ Efa) rijSiSj)(^, Sif- 

Here, we perform the same local gauge transformation. Then, we obtain 

A-1 V c::p(i^V. ^^ ^(/^,M)E.,.xiOxp(A:E,.,r,.-^.-v)(z;i,^V,)- 



1 ^ ^(%ill)<(f 5,.,)- (A6) 



Thus, from equations (A4) and (A6), we conclude that both sides of equation (6) 
coincide with each other. 



References 

[1] Bhatt R N and Young A P 1985 Phys. Rev. Lett. 54 924 

[2] Ogielski A T and Morgenstern I 1985 Phys. Rev. Lett. 54 928 

[3] Ogielski A T 1985 Phys. Rev. B 32 7384 

[4] Kawashima N and Young A P 1996 Phys. Rev. B 53 R484 

[5] Singh R R P and Chakravarty S 1986 Phys. Rev. Lett. 57 245 

[6] Bray A J and Moore M A 1984 J. Phys. C 17 L463 

[7] Swendsen R H and Wang J -S 1986 Phys. Rev. Lett. 57 2607 

[8] Bhatt R N and Young A P 1988 Phys. Rev. B 37 5606 

[9] Morgenstern I and Binder K 1980 Phys. Rev. B 22 288 

[10] Ozeki Y 1990 J. Phys. Soc. Jpn. 59 3531 

[11] Kawashima N and Ricgcr H 1997 Europhys. Lett. 39 85 

[12] Shirakura T and Matsubara F 1996 J. Phys. Soc. Jpn. 65 3138 

[13] Shiomi M, Matsubara F and Sliirakura T private communication 

[14] Shirakura T and Matsubara F 1997 Phys. Rev. Lett. 79 2887 

[15] Matsubara F, Shirakura T and Shiomi M 1998 Phys. Rev. B 58 R11821 

[16] Ozeki Y and Nishimori H 1987 J. Phys. Soc. Jpn. 56 3265 

[17] Kitatani H and Oguchi T 1992 J. Phys. Soc. Jpn. 61 1598 

[18] Nishimori H 1981 Prog. Theor. Phys. 66 1169 

[19] Kitatani H 1992 ,/. Phys. Soc. Jpn. 61 4049 

[20] Kitatani H and Seike T in preparation 

[21] Binder K 1991 Phys. Rev. Lett. 47 119 

[22] Kikuchi M and Ito N 1993 J. Phys. Soc. Jpn. 62 3052 



9 



300 



250 - 



Xs 

















L 




° 


□ 












_ 






□ 








X 7 




■ 


■ 










' 9 


- 






■ 


■ 


□ 




' 11 












■ 




♦ 13 






• • » 


• • 


• 






« 15 












» 




• 17 










* 






o 19 






is & 








* 


■ 21 






* * 




* 




* 


o 23 




















0.1 


0.2 


0.3 




0.4 


0.5 


0.6 



Figure 1. A plot of XSG versus 
T. 




0.9 ■ 



0.8 - 



0.7 - 



0.6 



r,=o.3 


J 


' °- „ v = 1.3 


" 7 


^* 77 = 0.225 


» 9 

« 11 




• 13 


\ 


» 15 


\ 


. 17 


* 


o 19 




■ 21 




□ 23 


• 









-1.5 



-0.5 0.5 1 1.5 2 2.5 



l/v 



Figure 4. A scaling plot for 
XSG- 



0.1 0.2 0.3 0.4 0.5 



0.75 



Figure 2. A plot of versus 
T. 





= 0.3 




I, 




i 


1 


V = 1.3 


X 7 
. 9 






« 11 




\ 


• 13 






• 15 






• 17 






• 19 






■ 21 






= 23 




i 






i 






5 





-1.5 -1 -Q.5 D.5 



1.5 2 2.5 



{T-TJL'" 
Figure 3. A scaling plot for gi^. 





% = 

v=2.6 

i I 


i 

X 7 

- 9 
» 11 
• 13 




1 

« s 


• 15 

• 17 






o 19 






■ 21 






= 23 




5 





Figure 5. A scaling plot for gj^, 
assuming Tc = 0. 
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Figure 6. A scaling plot for 
XSG> assuming Tc = 0. 



Figure 8. A scaling plot for 
Xsc; including the corrections to 
scaling with u)' = 0.5 and 6 = 
—0.3, assuming Tc = 0. 
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